#DESCRIPTION:
#Analysis of slope limits from results of factor of safety analysis


#1. libraries and modules

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import beta,gamma,gaussian_kde, kstest, ks_2samp, ttest_ind,norm
from lmfit import minimize, Parameters,report_errors
from scipy.interpolate import UnivariateSpline



#2. functions
def betadist(betaparams,B,pcF):
    #defining beta distribution parameters
    a=float(betaparams['a'].value)
    b=float(betaparams['b'].value)
    loc=float(betaparams['loc'].value)
    scale=float(betaparams['scale'].value)
    #creating fitted data
    model_pcF=beta.cdf(B,a,b,loc=loc,scale=scale)         
    #returning residual
    return (model_pcF-pcF)
    
def betadist_leastsquare_fitting(B,pcF,ai,bi,loc,scale):
    #defining power law variables
    betaparams = Parameters()
    betaparams.add('a', value=ai, vary=True, min=0.001, max=None)
    betaparams.add('b', value=bi, vary=True, min=0.001, max=None)
    betaparams.add('loc', value=loc, vary=True, min=0, max=loc)
    betaparams.add('scale', value=scale, vary=True, min=scale, max=5)
    #minimizing residuals
    result = minimize(betadist, betaparams, args=(B,pcF))
    print "RMSE: ", np.sqrt(np.mean(result.residual**2))
    Bfit = np.linspace(betaparams['loc'].value,betaparams['scale']+betaparams['loc'].value,90*10+1)
    pcFfit=beta.cdf(Bfit,
                    betaparams['a'].value,
                    betaparams['b'].value,
                    loc=betaparams['loc'].value,
                    scale=betaparams['scale'].value)
    
    return Bfit,pcFfit, np.sqrt(np.mean((result.residual)**2)),betaparams['a'].value,betaparams['b'].value,betaparams['loc'].value,betaparams['scale'].value


def power(powerparams,hl,B):
    hl=np.array(hl)
    B=np.array(B)
    #defining power distribution parameters
    a=float(powerparams['a'].value)
    b=float(powerparams['b'].value)
    #creating fitted data
    model_B=a*(hl**b)
    #returning residual
    return np.array((model_B-B))
    
def power_leastsquare_fitting(hl,B,ai,bi):
    #defining power law variables
    powerparams = Parameters()
    powerparams.add('a', value=ai, vary=True, min=0.001, max=None)
    powerparams.add('b', value=bi, vary=True, min=None, max=-0.00001)
    #minimizing residuals
    result = minimize(power, powerparams, args=(hl,B))
    hlfit = np.linspace(30,110,100)
    Bfit=powerparams['a'].value*(hlfit**powerparams['b'].value)
    
    return hlfit,Bfit, np.sqrt(np.mean((result.residual)**2))


def kde_minmode(data,x,max_num_mode,min_mode_pdf):
    kde=gaussian_kde(data)
    f=kde.factor
    f_list=np.linspace(f,(data.max()-data.min()),100)
    s=UnivariateSpline(x,kde(x),s=0)
    s1=UnivariateSpline(x,s(x,1),s=0)
    s2=UnivariateSpline(x,s1(x,1),s=0)
    extrema=s1.roots()
    
    maxima=extrema[np.where((s2(extrema)<0)*(s(extrema)>=min_mode_pdf))]
    
    if len(maxima)>max_num_mode:
        for q in range(1,len(f_list)):
            f=f_list[q]
            kde2=gaussian_kde(data,bw_method=f)
            s=UnivariateSpline(x,kde2(x),s=0)
            s1=UnivariateSpline(x,s(x,1),s=0)
            s2=UnivariateSpline(x,s1(x,1),s=0)
            extrema=s1.roots()
            maxima=extrema[np.where((s2(extrema)<0)*(s(extrema)>=min_mode_pdf))]
            if len(maxima)<=max_num_mode:
##                print 'modes: ',maxima
                break
        kde=gaussian_kde(data,bw_method=f)
##    else:
##        print maxima

    return kde,maxima


def plot_results(B,pcF,a,b,loc,scale,numsamp,sampsize,curcolor,lstyle,ax):
    Bfit,pcFfit,RMSE,a,b,loc,scale=betadist_leastsquare_fitting(B,pcF,a,b, loc, scale)
    beta_mode=scale*(a-1.)/(a+b-2.)

    ax[0].plot(Bfit,pcFfit,ls=lstyle,color=curcolor, label=Lith_list[j])
    ax[0].plot(B,pcF,'x',c=curcolor)
    ax[0].legend(fontsize='x-small', loc='best')
    if Lith_list[j]=='UM':
        ax[0].grid(b=None, which='major', axis='x')

    
    #generate random variables from beta distribution
    betarvs=beta.rvs(a,b,loc=loc,scale=scale,size=numsamp*100)
    ax[1].hist(betarvs, bins=np.linspace(0,5,1+int(5/0.05)),normed=True, histtype='step',lw=0.5,color=curcolor)
    ax[1].plot(Bfit, beta.pdf(Bfit,a,b,loc=loc,scale=scale),ls=lstyle,color=curcolor,lw=2)
    if Lith_list[j]=='UM':
        ax[1].grid(b=None, which='major', axis='x')
    

    samp_modes=np.empty(numsamp)
    for i in range(numsamp):
        indx=np.arange(len(betarvs))
        np.random.shuffle(indx)
        indx=indx[:sampsize]
        
        #computing and appending sample mode
        kde,maxima=kde_minmode(betarvs[indx],np.linspace(0,5,1+int(5/0.001)),2,0.1)
        samp_modes[i]=maxima[np.argmax(kde(maxima))]
        #ax[2].plot(np.linspace(0,5,1+int(5/0.001)),kde(np.linspace(0,5,1+int(5/0.001))),'-',c='0.6')

    
    kde,maxima=kde_minmode(samp_modes,np.linspace(0,5,1+int(5/0.001)),2,0.1)
    #ax[2].plot(np.linspace(0,5,1+int(5/0.001)),kde(np.linspace(0,5,1+int(5/0.001))),'-',c=curcolor)
    ksnorm_pval=kstest(samp_modes,'norm',args=(np.mean(samp_modes),np.std(samp_modes)))[1]
    if round(ksnorm_pval,3)>0.05:
        ax[2].hist(samp_modes, bins=np.linspace(0,5,1+int(5/0.05)),normed=True, histtype='step',lw=0.5,color=curcolor)
        ax[2].plot(np.linspace(0,5,1+int(5/0.001)),norm.pdf(np.linspace(0,5,1+int(5/0.001)),loc=np.mean(samp_modes),scale=np.std(samp_modes)),
                   '-',c=curcolor, lw=2)
        if Lith_list[j]=='UM':
            ax[2].grid(b=None, which='major', axis='x')
    
        

    return a,b, round(beta_mode,3),round(ksnorm_pval,3), samp_modes
    

def derivative(x,y,order):
    s=UnivariateSpline(x,y,s=0)
    dydx=s(x,order)
    return dydx
    
    

#3. MAIN
FSdata=pd.read_csv('R55_results_summary.csv',sep=',',header=0,index_col=None)

np.random.seed(12345)






hl_list=[110]#,50,70,90,110]
Lith_list=['CC','FC','LS', 'UM']
color_list=['b','g','r','c','m']

        
fig,ax=plt.subplots(nrows=3, ncols=1,sharex=True,sharey=False, figsize=(4.5,8))
main_plot=fig.add_subplot(111, frameon=False)
main_plot.set_xlabel('\nslope gradient, m/m',fontsize='small')
main_plot.set_xticklabels([])
main_plot.set_yticklabels([])
main_plot.set_xticks([])
main_plot.set_yticks([])

ax[0].set_ylabel('cumulative failure rate',fontsize='small')
plt.setp(ax[0].get_xticklabels(),fontsize='x-small')
plt.setp(ax[0].get_yticklabels(),fontsize='x-small')
ax[0].set_ylim(bottom=-0.02,top=1.02)

ax[1].set_ylabel('incremental failure rate',fontsize='small')
plt.setp(ax[1].get_xticklabels(),fontsize='x-small')
plt.setp(ax[1].get_yticklabels(),fontsize='x-small')

ax[2].set_ylabel('mode pdf',fontsize='small')
plt.setp(ax[2].get_xticklabels(),fontsize='x-small')
plt.setp(ax[2].get_yticklabels(),fontsize='x-small')

df1=pd.DataFrame()
df2=pd.DataFrame()

samp_modes_list=[]
for i in range(len(hl_list)):

    
    for j in range(len(Lith_list)):
               
        cur_dat=FSdata[(FSdata.hl==hl_list[i])&(FSdata.Lith==Lith_list[j])]
        #convert from degrees to m/m
        B=np.tan(np.radians((cur_dat['B'].values)[::-1]))
        pcF=(cur_dat['pc_failure'].values)[::-1]/100.
        toe_exit=(cur_dat['toe exit'].values)[::-1]
                
        curcolor=color_list[j]

        #initializing Beta distribution parameters for fitting
        a,b,loc,scale=5,2,0.00001,4.99999
        
        # Setting pcF=100% for slope face exits (toe exit>10)
        pcF2=np.where(toe_exit<10,pcF,np.ones(len(pcF)))
        dfrow={'bl':hl_list[i], 'Li':Lith_list[j],
               '15':np.round(pcF2[0]*100,1),
               '30':np.round(pcF2[1]*100,1),
               '45':np.round(pcF2[2]*100,1),
               '60':np.round(pcF2[3]*100,1),
               '75':np.round(pcF2[4]*100,1)}
        df1=df1.append(dfrow, ignore_index=True)

        numsamp=30
        sampsize=30
        a,b, beta_mode,ksnorm_pval,samp_modes=plot_results(B,pcF2,a,b,loc,scale,numsamp,sampsize,curcolor,'-',ax)
        samp_modes_list.append(samp_modes)
        dfrow={'Li':Lith_list[j],
               'Beta_mode':beta_mode,'Beta_a':a, 'Beta_b':b,
               'Beta_loc':round(loc,3), 'Beta_scale':round(scale,3),
               'SDmode_N':numsamp,'SDmode_ksnorm_pval':ksnorm_pval, 'SDmode_mean':np.mean(samp_modes), 'SDmode_std':np.std(samp_modes)}
        df2=df2.append(dfrow, ignore_index=True)

##        if i==0:
##            curax.set_title(Lith_list[j],fontsize='small')
       

            



df1.set_index(['Li','bl'],inplace=True,drop=True)
df2.set_index(['Li'],inplace=True,drop=True)
print df1
print df2
#df1.to_csv('r55_df1_pcF_analysis_results.csv',sep=',')
#df2.to_csv('r55_df2_pcF_analysis_results.csv',sep=',')

fig.tight_layout()

df3=pd.DataFrame()
for k in range(len(samp_modes_list)):
    for l in range(k+1,len(samp_modes_list)):
        cur=samp_modes_list[k]
        comp=samp_modes_list[l]

        t,pval=ttest_ind(cur,comp,equal_var=False)
        ks,kspval=ks_2samp(cur,comp)

        dfrow={'li_cur':Lith_list[k],
               'li_comp':Lith_list[l],
               't_pval': round(pval,3),
               'diff_mean':round(np.mean(cur)-np.mean(comp),3),
               'diff_median':round(np.median(cur)-np.median(comp),3),
               'diff_std':round(np.std(cur)-np.std(comp),3),
               'kspval':round(kspval,3)}
        df3=df3.append(dfrow, ignore_index=True)

df3.set_index(['li_cur','li_comp'],inplace=True,drop=True)
#df3.to_csv('r55_df3_pcF_analysis_results.csv',sep=',')
print df3
#plt.savefig('r55_pcF_analysis_results.png')
        



##
##fig,ax=plt.subplots(nrows=2,ncols=1, sharex=True, sharey=True, figsize=(3.5,6))
##main_plot=fig.add_subplot(111, frameon=False)
##main_plot.set_ylabel('slope gradient, degrees\n\n ')
##main_plot.set_xlabel('\n\nbase length, m')
##main_plot.set_xticklabels([])
##main_plot.set_yticklabels([])
##main_plot.set_xticks([])
##main_plot.set_yticks([])
##for j in range(len(Lith_list)):
##   
##    
##    cur_df2=df2[(df2.Li==Lith_list[j])]
##    hl=cur_df2.hl.astype(float).values
##    dydx=cur_df2.f1_max.astype(float).values
##    d2ydx2=cur_df2.f2_max.astype(float).values
##
##    #plt.plot(hl, dydx, 'ko')
##    #plt.plot(hl,d2ydx2, 'kx')
##
##    plt.sca(ax[0])
##
##    ai,bi=30,-1
##
##    line,=plt.plot(hl,dydx, '.')
##    curcolor=plt.getp(line, 'color')
##    
##    
####    hl2=np.linspace(hl.min(),hl.max(),100)
####    s1=UnivariateSpline(hl,dydx,s=1,k=2)
####    dydx=s1(hl2)
####    plt.plot(hl2,dydx, ':', color=curcolor)
##
##    hl2,dydx,RMSE=power_leastsquare_fitting(hl,cur_df2.f1_max.astype(float).values,ai,bi)
##    plt.plot(hl2, dydx, '-', color=curcolor, label=Lith_list[j])
##    
##    if j==len(Lith_list)-1:
##        plt.title('maximum pcF gradient', fontsize='small')
##        plt.legend(fontsize='x-small', loc='upper right')
##        ax[0].set_yticks((15,30,45,60,75))
##
##    plt.sca(ax[1])
##    plt.plot(hl,d2ydx2, '.', color=curcolor)
##
####    s2=UnivariateSpline(hl,d2ydx2, s=1,k=2)
####    d2ydx2=s2(hl2)
####    plt.plot(hl2,d2ydx2, ':', color=curcolor)
##
##    hl2,d2ydx2,RMSE=power_leastsquare_fitting(hl,cur_df2.f2_max.astype(float).values,ai,bi)
##    plt.plot(hl2,d2ydx2, '-', color=curcolor)
##
##    if j==len(Lith_list)-1:
##        plt.title('maximum pcF curvature', fontsize='small')
##        ax[1].set_xticks((30,50,70,90,110))
##        ax[1].set_yticks((15,30,45,60,75))
##        ax[1].set_ylim(15,75)
##        ax[1].set_xlim(25,115)
##
##fig.tight_layout()


##        
##fig,ax=plt.subplots(nrows=2,ncols=2, sharex=True, sharey=True, figsize=(5,6))
##main_plot=fig.add_subplot(111, frameon=False)
##main_plot.set_ylabel('slope gradient, degrees\n\n  ')
##main_plot.set_xlabel('\n\nbase length, m')
##main_plot.set_xticklabels([])
##main_plot.set_yticklabels([])
##main_plot.set_xticks([])
##main_plot.set_yticks([])
##j2=[(0,0),(0,1),(1,0),(1,1)]
##for j in range(len(Lith_list)):
##   
##    
##    cur_df2=df2[(df2.Li==Lith_list[j])]
##    hl=cur_df2.hl.astype(float).values
##
##    pc10=cur_df2['10'].astype(float).values
##    pc25=cur_df2['25'].astype(float).values
##    pc50=cur_df2['50'].astype(float).values
##    pc75=cur_df2['75'].astype(float).values
##    pc100=cur_df2['100'].astype(float).values
##    
####    hl2=np.linspace(hl.min(),hl.max(),100)
####    s10=UnivariateSpline(hl,pc10,s=1,k=2)
####    s25=UnivariateSpline(hl,pc25,s=1,k=2)
####    s50=UnivariateSpline(hl,pc50,s=1,k=2)
####    s75=UnivariateSpline(hl,pc75,s=1,k=2)
####    s100=UnivariateSpline(hl,pc100,s=1,k=2)
##
##    hl2,s10,RMSE=power_leastsquare_fitting(hl,pc10,ai,bi)
##    hl2,s25,RMSE=power_leastsquare_fitting(hl,pc25,ai,bi)
##    hl2,s50,RMSE=power_leastsquare_fitting(hl,pc50,ai,bi)
##    hl2,s75,RMSE=power_leastsquare_fitting(hl,pc75,ai,bi)
##    hl2,s100,RMSE=power_leastsquare_fitting(hl,pc100,ai,bi)
##
##   
##    curax=ax[0,0]
##    if j==len(Lith_list)-1:
##        curax.set_yticks((15,30,45,60,75,90))
##        curax.set_xticks((30,50,70,90,110))
##        curax.set_ylim(15,90)
##        curax.set_xlim(25,115)
##    line,=curax.plot(hl2, s25, '-', label=Lith_list[j])
##    curcolor=plt.getp(line,'color')
##    curax.plot(hl, pc25, '.', color=curcolor)
##    curax.set_title('25% failure',fontsize='small')
##    curax.legend(loc='upper right',fontsize='x-small')
##
##
##    curax=ax[0,1]
##    if j==len(Lith_list)-1:
##        curax.set_yticks((15,30,45,60,75,90))
##        curax.set_xticks((30,50,70,90,110))
##        curax.set_ylim(15,90)
##        curax.set_xlim(25,115)
##    line,=curax.plot(hl2, s50, '-', label=Lith_list[j])
##    curcolor=plt.getp(line,'color')
##    curax.plot(hl, pc50, '.', color=curcolor)
##    curax.set_title('50% failure',fontsize='small')
##
##    curax=ax[1,0]
##    if j==len(Lith_list)-1:
##        curax.set_yticks((15,30,45,60,75,90))
##        curax.set_xticks((30,50,70,90,110))
##        curax.set_ylim(15,90)
##        curax.set_xlim(25,115)
##    line,=curax.plot(hl2, s75, '-', label=Lith_list[j])
##    curcolor=plt.getp(line,'color')
##    curax.plot(hl, pc75, '.', color=curcolor)
##    curax.set_title('75% failure',fontsize='small')
##
##
##    curax=ax[1,1]
##    if j==len(Lith_list)-1:
##        curax.set_yticks((15,30,45,60,75,90))
##        curax.set_xticks((30,50,70,90,110))
##        curax.set_ylim(15,90)
##        curax.set_xlim(25,115)
##    line,=curax.plot(hl2, s100, '-', label=Lith_list[j])
##    curcolor=plt.getp(line,'color')
##    curax.plot(hl, pc100, '.', color=curcolor)
##    curax.set_title('100% failure',fontsize='small')  
##
##
##        
##fig.tight_layout()
##print df2
plt.show()







    
